##########################################################
# Código para a análise multinomial
# 
# Autoria: Leonardo Rodrigues
# Código escrito e executado no RStudio (Version 1.3.1073)
# Tempo aproximado para rodar o script: 3.3 horas
##########################################################
start_time <- Sys.time()

### 1 Preparação ----
# Leitura dos pacotes

library(webshot)
library(dplyr)
library(stringr)
library(tidyr)
library(tools)
library(summarytools)
library(stargazer)
library(mice)
library(miceadds)
library(nnet)


#Leitura do banco
enade <- read.table("data/originaisEnade/enade.txt", sep="\t")

#Mudando nome das categorias para aproveitar script anterior
#estabelecendo a categoria de referência para o modelo

enade$cor = ifelse(enade$cor == "Brancos", "Branca",
                   ifelse(enade$cor == "Negros", "Não-Branca", NA))

enade$year <- as.factor(enade$ciclo)
enade$year = relevel(enade$year, ref = "2")
enade$tese = ifelse(str_detect(enade$curso,"ENGENHARIA"), enade$curso, "outros")
enade11 <- enade
enade11$year <- enade11$ciclo
enade11$year <- as.factor(enade11$year)
enade11$year = relevel(enade11$year, ref = "2")
enade11$grupo2 <- enade11$diploma
enade11$grupo2 <- as.factor(enade11$grupo2)
enade11$grupo2 = relevel(enade11$grupo2, ref = "Engenharias")
enade11$regiao <- enade11$regiao1
#excluir banco que não será utilizado
enade <- NULL 

#selecionando variáveis de interesse
enade11 <- enade11[, c("grupo2", "year", "age", "regiao", "sexo", "cor", "edu")]

#transformando variáveis em fator
enade11$edu <- as.factor(enade11$edu)
enade11$cor <- as.factor(enade11$cor)
enade11$sexo <- as.factor(enade11$sexo)
enade11$regiao <- as.factor(enade11$regiao)
enade11$age <- as.factor(enade11$age)

### 2 Imputação de valores para os missing ----
#Separando os bancos
enade1 <- enade11[enade11$year==2,]
enade2 <- enade11[enade11$year==4,]
rm(enade11)

#fazendo a imputação (uma imputação para cada ciclo)
#Imputação ciclo 1
imputed_data1 <- mice(enade1, m=5, maxit=25, method = 'logreg')

#Imputação ciclo 2
imputed_data2 <- mice(enade2, m=5, maxit=25, method = 'logreg')

#juntando os bancos de dados imputados
imputed_total <- rbind(imputed_data1, imputed_data2)

### 3 Modelo Multinomial ----
max_iterations <- 500

multil <- with(imputed_total, multinom(grupo2 ~ age + regiao + sexo*year + cor*year + edu*year, maxit = max_iterations))

multilcomb <- pool(multil) #juntando os cinco modelos

output <- summary(multilcomb) #resumo da junção

#Transformar os resultados em um objeto "MULTINOM"
# Copiando um dos modelos para utilizar a estrutura semelhante ao objeto "Multinom"
pooled = multil$analyses[[1]]

#padronizar o modelo "Pooled" com o objeto R "Multinom"
teste1  <- rep(0, 12) #criando um vetor com 0. O modelo Multinom inclui o 0 nos coeficientes da categoria padrão (engenharias)
teste2  <- output$estimate #colando os coeficientes do modelo Pooled em um novo vetor

vetor <- c(teste1, teste2) #juntando os vetores no mesmo padrão do objeto "pooled".
vetor <- append(vetor, 0, 22)
vetor <- append(vetor, 0, 33)
vetor <- append(vetor, 0, 44)
vetor <- append(vetor, 0, 55)

# Substituir os coeficientes ajustados pelo "pooled estimates"
pooled$wts <- vetor #colando os coeficientes do modelo Pooled em um objeto "Multinom"
#as probabilidades preditas serão feitas a partir dos coeficientes do modelo "Pooled"

multil.rrr = exp(coef(pooled))
multil <- pooled

head(pp <- fitted(multil))

#Obter a razão entre as propabilidades preditas (para todas as áreas)
#criando um data.frame para salvar a razão das probabilidades preditas ("PPS")
pps <- data.frame("area" = c("Engenharias", "Bacharelados", "Direito", "Licenciaturas", "Medicina", "Tecnólogos"),
                  "MF1" = NA, "FM1" = NA, "SM1" = NA, "MS1" = NA, "BN1" = NA, "NB1" = NA,
                  "MF2" = NA, "FM2" = NA, "SM2" = NA, "MS2" = NA, "BN2" = NA, "NB2" = NA)

#criando um data.frame para salvar as probabilidades preditas relacionadas a escolaridade dos pais, mantendo as outras variáveis em suas médias
mediaEDU <- data.frame(edu = c("Ensino superior ou mais","Até Ensino Médio", "Ensino superior ou mais", "Até Ensino Médio"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Feminino", "Feminino", "Feminino", "Feminino"),
                       cor = c("Branca","Branca", "Branca","Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

#obtendo as probabilidades preditas para a variável escolaridade dos pais
mediaEDU[, c("pred.prob")] <- predict(multil, newdata=mediaEDU, type="probs")

#obtendo a razão das probabilidades preditas entre as categorias da variável escolaridade dos pais para cada ciclo.
ratEDU1 <- t(data.frame(mediaEDU[1,7]/mediaEDU[2,7])) #SM1
ratEDU2 <- t(data.frame(mediaEDU[2,7]/mediaEDU[1,7])) #MS1
ratEDU3 <- t(data.frame(mediaEDU[3,7]/mediaEDU[4,7])) #SM2
ratEDU4 <- t(data.frame(mediaEDU[4,7]/mediaEDU[3,7])) #MS2

#salvando a razão das probabilidades preditas no data.frame "PPS"
pps[,"SM1"] <- ratEDU1
pps[,"MS1"] <- ratEDU2
pps[,"SM2"] <- ratEDU3
pps[,"MS2"] <- ratEDU4

#O mesmo processo foi utilizado para a variável de Sexo
mediaSEXO <- data.frame(edu = c("Até Ensino Médio","Até Ensino Médio", "Até Ensino Médio", "Até Ensino Médio"),
                        year = c("2", "2", "4", "4"),
                        sexo = c("Masculino", "Feminino", "Masculino", "Feminino"),
                        cor = c("Branca","Branca", "Branca","Branca"),
                        age = c("20-30", "20-30", "20-30", "20-30"),
                        regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaSEXO[, c("pred.prob")] <- predict(multil, newdata=mediaSEXO, type="probs")

ratSEXO1 <- t(data.frame(mediaSEXO[1,7]/mediaSEXO[2,7])) #MF1
ratSEXO2 <- t(data.frame(mediaSEXO[2,7]/mediaSEXO[1,7])) #FM1
ratSEXO3 <- t(data.frame(mediaSEXO[3,7]/mediaSEXO[4,7])) #MF2
ratSEXO4 <- t(data.frame(mediaSEXO[4,7]/mediaSEXO[3,7])) #FM2

pps[,"MF1"] <- ratSEXO1
pps[,"FM1"] <- ratSEXO2
pps[,"MF2"] <- ratSEXO3
pps[,"FM2"] <- ratSEXO4

#O mesmo processo foi utilizado para a variável de Cor/Raça
mediaCOR <- data.frame(edu = c("Até Ensino Médio","Até Ensino Médio", "Até Ensino Médio", "Até Ensino Médio"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Feminino", "Feminino", "Feminino", "Feminino"),
                       cor = c("Branca","Não-Branca", "Branca","Não-Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaCOR[, c("pred.prob")] <- predict(multil, newdata=mediaCOR, type="probs")

ratCOR1 <- t(data.frame(mediaCOR[1,7]/mediaCOR[2,7])) #BN1
ratCOR2 <- t(data.frame(mediaCOR[2,7]/mediaCOR[1,7])) #NB1
ratCOR3 <- t(data.frame(mediaCOR[3,7]/mediaCOR[4,7])) #BN2
ratCOR4 <- t(data.frame(mediaCOR[4,7]/mediaCOR[3,7])) #NB2

pps[,"BN1"] <- ratCOR1
pps[,"NB1"] <- ratCOR2
pps[,"BN2"] <- ratCOR3
pps[,"NB2"] <- ratCOR4

#Obtendo e salvando a razão entre as probabilidades preditas em que a referência é a categoria "dominante".
#Por exemplo, se a chance de homens serem concluintes nas licenciaturas, em comparação com as outras áreas, é menor que a de mulheres
#o resultado foi multiplicado por -1, para indicar a direção da desigualdade.
pps$sexo1 <- ifelse(pps$MF1 > 1, pps$MF1,
                    ifelse(pps$MF1 < 1, (pps$FM1)*-1, NA))
pps$sexo2 <- ifelse(pps$MF2 > 1, pps$MF2,
                    ifelse(pps$MF2 < 1, (pps$FM2)*-1, NA))

pps$edu1 <- ifelse(pps$SM1 > 1, pps$SM1,
                   ifelse(pps$SM1 < 1, (pps$MS1)*-1, NA))
pps$edu2 <- ifelse(pps$SM2 > 1, pps$SM2,
                   ifelse(pps$SM2 < 1, (pps$MS2)*-1, NA))

pps$cor1 <- ifelse(pps$BN1 > 1, pps$BN1,
                   ifelse(pps$BN1 < 1, (pps$NB1)*-1, NA))
pps$cor2 <- ifelse(pps$BN2 > 1, pps$BN2,
                   ifelse(pps$BN2 < 1, (pps$NB2)*-1, NA))


#salvando os resultados para fazer os gráficos.
write.table(pps, "data/output/ratioppsareas.txt", sep="\t")

# Obtendo a razão entre as probabilidades preditas para cada variável (áreas imperiais)
#criando banco para o Ciclo 1
media1 <- data.frame(sexo = c("Masculino", "Feminino", "Masculino", "Masculino", "Feminino", "Feminino", "Masculino", "Feminino"),
                     edu = c("Ensino superior ou mais", "Ensino superior ou mais", "Até Ensino Médio", "Ensino superior ou mais", "Até Ensino Médio",
                             "Ensino superior ou mais", "Até Ensino Médio", "Até Ensino Médio"),
                     cor = c("Branca","Branca", "Branca","Não-Branca", "Branca", "Não-Branca", "Não-Branca", "Não-Branca"),
                     age = c("20-30", "20-30", "20-30", "20-30", "20-30", "20-30", "20-30", "20-30"),
                     regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste", "Sul e Sudeste", "Sul e Sudeste", "Sul e Sudeste"),
                     year = c("2", "2", "2", "2", "2", "2", "2", "2"))

#criando banco para o Ciclo 2
media2 <- data.frame(sexo = c("Masculino", "Feminino", "Masculino", "Masculino", "Feminino", "Feminino", "Masculino", "Feminino"),
                     edu = c("Ensino superior ou mais", "Ensino superior ou mais", "Até Ensino Médio", "Ensino superior ou mais", "Até Ensino Médio",
                             "Ensino superior ou mais", "Até Ensino Médio", "Até Ensino Médio"),
                     cor = c("Branca","Branca", "Branca","Não-Branca", "Branca", "Não-Branca", "Não-Branca", "Não-Branca"),
                     age = c("20-30", "20-30", "20-30", "20-30", "20-30", "20-30", "20-30", "20-30"),
                     regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste", "Sul e Sudeste", "Sul e Sudeste", "Sul e Sudeste"),
                     year = c("4", "4", "4", "4", "4", "4", "4", "4"))


#obtendo as probabilidades preditas para todas as combinações das variáveis ciclo 1
media1[, c("pred.prob")] <- predict(multil, newdata=media1, type="probs")

#obtendo as probabilidades preditas para todas as combinações das variáveis ciclo 2
media2[, c("pred.prob")] <- predict(multil, newdata=media2, type="probs")

#backup da pps
write.table(media1, "data/output/ppsinternaciclo1.txt", sep="\t")
write.table(media2, "data/output/ppsinternaciclo2.txt", sep="\t")

#Retirando as variáveis que não serão utilizadas (região e idade)
media1[4:5] <- list(NULL)
media2[4:5] <- list(NULL)

#selecionando apenas as variáveis relacionadas às probabilidades preditas
pp1 <- as.data.frame(media1$pred.prob)
pp2 <- as.data.frame(media2$pred.prob)

#Fazendo a razão das probabilidades preditas do perfil dominante pelos demais perfis
pp1 <- pp1 %>% 
  mutate(ppEngenharias = pp1[1,1]/Engenharias,
         ppDireito = pp1[1,3]/Direito,
         ppMedicina = pp1[1,5]/Medicina) %>%
  select(ppEngenharias, ppDireito, ppMedicina)

pp2 <- pp2 %>% 
  mutate(ppEngenharias = pp2[1,1]/Engenharias,
         ppDireito = pp2[1,3]/Direito,
         ppMedicina = pp2[1,5]/Medicina) %>%
  select(ppEngenharias, ppDireito, ppMedicina)

#montando um data frame com a razão das probabilidades preditas do perfil dominante pelos demais perfis
media1 <- media1 %>% unite(comb, c("sexo", "edu", "cor"), sep = " ") %>%
  select(comb, year)

media2 <- media2 %>% unite(comb, c("sexo", "edu", "cor"), sep = " ") %>%
  select(comb, year)

pp1 <- merge(media1, pp1, by=0)
pp2 <- merge(media2, pp2, by=0)
pps <- rbind(pp1, pp2)

pps$comb <- dplyr::recode(pps$comb,"Masculino Ensino superior ou mais Branca" = "Homem Ensino Superior Branco",
                          "Feminino Ensino superior ou mais Branca" = "Mulher Ensino Superior Branca",
                          "Masculino Até Ensino Médio Branca" = "Homem Ensino Médio Branco",
                          "Masculino Ensino superior ou mais Não-Branca" = "Homem Ensino Superior Negro",
                          "Feminino Até Ensino Médio Branca" = "Mulher Ensino Médio Branca",
                          "Feminino Ensino superior ou mais Não-Branca" = "Mulher Ensino Superior Negra",
                          "Masculino Até Ensino Médio Não-Branca" = "Homem Ensino Médio Negro",
                          "Feminino Até Ensino Médio Não-Branca" = "Mulher Ensino Médio Negra")

pps$year <- dplyr::recode(pps$year, "2" = "1",
                          "4" = "2")

#salvando o banco de dados para fazer os gráficos 
write.table(pps, "data/output/ppsinterna.txt", sep="\t")

#End of script
end_time <- Sys.time()
end_time - start_time #3.320241 hours